Fitting models to surface

In his series of papers studying and comparing various magnetopause empirical models, [Nguyen et al., 2022] defines numerically these models in a publicly available script. We will use some of these functions to understand whether they could describe the simulation’s magnetopause. We will be initializing each of the following fits with the simulation parameters as indicated in Initial OpenGGCM simulation parameters.. For the parameter estimation we will use directly the built-in curve-fit function of SciPy.

Here we try fitting the Shue model to the full 3D data of the maximum emissivity surface:

Show code cell source Hide code cell source

theta_grid, phi_grid = Theta[:,0,:], Phi[:,0,:]
gamma = 0
P_dyn = 12.5
P_m = 1.5
B_z = 5
Omega = np.radians(0)

phi_flat = np.ravel(phi_grid)
theta_flat = np.ravel(theta_grid)
r_flat = np.ravel(r_surface)

# Mask to keep only finite values
mask = np.isfinite(theta_flat) & np.isfinite(r_flat)
theta_flat_valid = theta_flat[mask]
r_flat_valid = r_flat[mask]
popt, pcov = curve_fit(model.MP_Shue1998, theta_flat_valid, r_flat_valid, p0=[12.5, 5])

perr = np.sqrt(np.diag(pcov))
param_names = ['Pd', 'Bz']
for name, value, error in zip(param_names, popt, perr):
    print(f"{name} = {value:.4f} ± {error:.4f}")

DP = popt[0]
Bz = popt[1]
r0 = (10.22+1.29*np.tanh(0.184*(Bz+8.14)))*DP**(-1./6.6)
a = (0.58-0.007*Bz)*(1+0.024*np.log(DP))
print(f"r0={r0}, a={a}")

from sklearn.metrics import mean_squared_error

# Compute Shue model surface with fitted parameters
r_shue_fit = model.MP_Shue1998(theta_grid, popt[0], popt[1])

# Compute RMSE (mask out invalid values if needed)
mask = ~np.isnan(r_surface) & ~np.isnan(r_shue_fit)
rmse = np.sqrt(mean_squared_error(r_surface[mask], r_shue_fit[mask]))
print(f"RMSE between r_surface and fitted Shue surface: {rmse} RE")

Pd = 11.3580 ± 0.0127
Bz = 6.5704 ± 0.1762
r0=7.956991679798463, a=0.5651498286559324
RMSE between r_surface and fitted Shue surface: 0.5972698393048848 RE

Show code cell source Hide code cell source

theta = np.arange(0, np.pi/2, 0.01)
phi = np.arange(0, 2*np.pi, 0.01)
P_dyn = popt[0]
B_z = popt[1]
r_shue = model.MP_Shue1998(theta_grid, P_dyn, B_z)
x, y, z = model.get_cartesian(r_shue, theta_grid, phi_grid)
fig = go.Figure(data=[go.Surface(
    x=x,
    y=y,
    z=z,
    surfacecolor=r_shue,
    colorscale='thermal',
    cmin=0,
    cmax=np.nanmax(r_shue),
    colorbar=dict(title="r value")
)])
fig.update_layout(
    title="Shue model",
    template="plotly_dark",
    scene=dict(
        xaxis_title="X",
        yaxis_title="Y",
        zaxis_title="Z"
    )
)
fig.show()

As expected, the Shue model cannot accurately describe the surface, rather it fits to the cusps instead of capturing the outer boundary.

The same process was implemented for the Lin model, which accounts for a significant indentation due to the cusps. However, as seen in \autoref{fig:lin_fit} the indentation is a localized feature that does not seem to describe the geometry of the extracted surface. The pressures do not converge, and we end up with a significant error both in their estimation and the final \ac{RMSE}.

Show code cell source Hide code cell source

P_dyn = 12.5
P_m = 0
B_z = -5

theta_fit, phi_fit = Theta[:,0,:], Phi[:,0,:]
phi_flat = np.ravel(phi_fit)
theta_flat = np.ravel(theta_fit)
r_flat = np.ravel(r_surface)



def wrapper(coords, Pd, Pm, Bz):
    phi, theta = coords
    return model.MP_Lin2010(phi, theta, Pd, Pm, Bz)
# Initial guess for the parameters
initial_guess = [P_dyn, P_m, B_z]
bounds = ([0,0,-10], [np.inf, np.inf, 10])  
# Fit the model
popt, pcov = curve_fit(wrapper, (phi_flat, theta_flat), r_flat, p0=initial_guess, bounds=bounds)

# popt contains the best-fit parameters
perr = np.sqrt(np.diag(pcov))
param_names = ['Pd', 'Pm', 'Bz']
for name, value, error in zip(param_names, popt, perr):
    print(f"{name} = {value:.4f} ± {error:.4f}")

Pd = popt[0]
Pm = popt[1]
Bz = popt[2]
tilt = 0 
Bx = 0 
By = 0
P = Pd+Pm

def exp2(i):
    return a[i]*(np.exp(a[i+1]*Bz)-1)/(np.exp(a[i+2]*Bz)+1)
a = [12.544,
         -0.194,
         0.305,0.0573,2.178]

r0 = a[0]*P**a[1]*(1+exp2(2))

# r0 = (10.56+0.956*np.tanh(0.1795*(Bz+10.78)))*P**(-0.1699)

print(f"r0={r0}, a={a}")

# Compute Shue model surface with fitted parameters

r_lin_fit = model.MP_Lin2010(phi_grid, theta_grid, popt[0], popt[1], popt[2])

# Compute RMSE (mask out invalid values if needed)
mask = ~np.isnan(r_surface) & ~np.isnan(r_lin_fit)
rmse = np.sqrt(mean_squared_error(r_surface[mask], r_lin_fit[mask]))
print(f"RMSE between r_surface and fitted Shue surface: {rmse}")

Pd = 4.0897 ± 18537.0867
Pm = 0.0003 ± 18537.0866
Bz = -7.8077 ± 0.0817
r0=8.494669436219956, a=[12.544, -0.194, 0.305, 0.0573, 2.178]
RMSE between r_surface and fitted Shue surface: 0.5835083284032755

Show code cell source Hide code cell source

theta = np.arange(0, np.pi/2, 0.01)
phi = np.arange(0, 2*np.pi, 0.01)
theta_full, phi_full = np.meshgrid(theta, phi, indexing='ij')
r = model.MP_Lin2010(
    th_in=theta_full,
    phi_in=phi_full,
    Pd = popt[0],
    Pm = popt[1],
    Bz = popt[2],
    tilt =0)

x, y, z = model.get_cartesian(r, theta_full, phi_full)

fig = go.Figure(data=[go.Surface(
    x=x,
    y=y,
    z=z,
    surfacecolor=r,
    colorscale='thermal',
    cmin=0,
    cmax=np.nanmax(r),
    colorbar=dict(title="r value")
)])

fig.update_layout(
    title="Lin model",
    template="plotly_dark",
    scene=dict(
        xaxis_title="X",
        yaxis_title="Y",
        zaxis_title="Z"
    )
)

fig.show()

Finally, we attempted to fit the Liu model, which is a more resembling geometry of the original surface. The model could not predict the positions of the cusps to be so close to the subsolar point, and ended up fitting the entire \([-65,65]\) degrees domain of our data as the front lobe of the magnetopause.

Show code cell source Hide code cell source

P_dyn = 12.5
B_z = -5
P_m = 0

def wrapper(coords, Pd, Pm, Bz):
    phi, theta = coords
    return model.MP_Liu2015(phi, theta, Pd, Pm, 0, 0, Bz, 0)
# Initial guess for the parameters

initial_guess = [P_dyn, P_m, B_z]
bounds = ([0,0,-10], [np.inf, np.inf, 10])  
# Fit the model
popt, pcov = curve_fit(wrapper, (phi_flat, theta_flat), r_flat, p0=initial_guess, bounds=bounds)

# popt contains the best-fit parameters
perr = np.sqrt(np.diag(pcov))
param_names = ['Pd', 'Pm', 'Bz']
for name, value, error in zip(param_names, popt, perr):
    print(f"{name} = {value:.4f} ± {error:.4f}")


r_liu_fit = model.MP_Liu2015(phi_grid, theta_grid, popt[0], popt[1], 0 , 0, popt[2],0)

Pd = popt[0]
Pm = popt[1]
Bz = popt[2]
tilt = 0 
Bx = 0 
By = 0
P = Pd+Pm

r0 = (10.56+0.956*np.tanh(0.1795*(Bz+10.78)))*P**(-0.1699)


print(f"r0={r0}, a={a}")

# Compute RMSE (mask out invalid values if needed)
mask = ~np.isnan(r_surface) & ~np.isnan(r_lin_fit)
rmse = np.sqrt(mean_squared_error(r_surface[mask], r_liu_fit[mask]))
print(f"RMSE between r_surface and fitted Shue surface: {rmse}")

Pd = 6.9549 ± 1.3755
Pm = 0.0000 ± 1.3737
Bz = -1.0846 ± 0.0426
r0=8.24208984586283, a=[12.544, -0.194, 0.305, 0.0573, 2.178]
RMSE between r_surface and fitted Shue surface: 0.6109231997180545

Show code cell source Hide code cell source

P_dyn = popt[0]
B_z = popt[2]
P_m = popt[1]

theta = np.arange(0, np.pi/2, 0.05)
phi = np.arange(0, 2*np.pi+0.05, 0.05)
theta_grid, phi_grid = np.meshgrid(theta, phi, indexing='ij')

r_liu = model.MP_Liu2015(
    th_in=theta_grid,
    phi_in=phi_grid,
    Pd = P_dyn,
    Pm = P_m,
    Bx = 0,
    By= 0,
    Bz = B_z,
    tilt =0)

x, y, z = model.get_cartesian(r_liu, theta_grid, phi_grid)

fig = go.Figure(data=[go.Surface(
    x=x,
    y=y,
    z=z,
    surfacecolor=r_liu,
    colorscale='thermal',
    cmin=0,
    cmax=np.nanmax(r),
    colorbar=dict(title="r value")
)])

fig.update_layout(
    title="Liu2015 model",
    template = 'plotly_dark',
    scene=dict(
        xaxis_title="X",
        yaxis_title="Y",
        zaxis_title="Z"
    )
)

fig.show()

Even if we ignore the shape differences of the models with our simulation data as an approximation, the large contribution from the jumps of the maximum emissivity does not allow for reliable fits. All the above models, when fitted, overestimated the subsolar distance, since the contribution of the jumps is non-negligible. Simultaneously, errors of the order of \(1\,RE\) cannot be neglected when we are trying to verify the tangent hypothesis, in order to utilize it for sub-RE detection of the magnetopause.